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Abstract 

We introduce an approach for the accurate calculation of thermal properties of classical 
nanoclusters. Based on a recently developed enhanced sampling technique, replica exchange 
metadynamics, the method yields the true free energy of each relevant cluster structure, di- 
rectly sampling its basin and measuring its occupancy in full equilibrium. All entropy sources, 
whether vibrational, rotational anharmonic and especially configurational - the latter often for- 
gotten in many cluster studies - are automatically included. For the present demonstration we 
choose the water nonamer (H20)9, an extremely simple cluster which nonetheless displays a 
sufficient complexity and interesting physics in its relevant structure spectrum. Within a stan- 
dard TIP4P potential description of water, we find that the nonamer second relevant structure 
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possesses a higher configurational entropy than the first, so that the two free energies surpris- 
ingly cross for increasing temperature. 

Introduction 

The physics of clusters formed by a small number of atoms or molecules is basic to many branches 
of science.^ Clusters are the primary nanosystem where matter behaves differently from its usual 
bulk aggregation. In condensed matter physics it has long been of interest how collective features 
of macroscopic systems, such as ordering phenomena and phase transitions, emerge in a nutshell 
in small size clusters.-^ As is well known, it often proves harder to predict the properties of a small 
aggregate than those of an infinitely large system, where one can take advantage of regularity. In 
computational physics, identifying and characterizing cluster structures provide a prototypically 
difficult optimization problem.^ The potential energy landscape as a function of the atomic or 
molecular coordinates is generally complex, with many local minima separated by large barriers 
even in relatively small clusters, where the simple identification of the lowest energy structure may 
be a serious computational challenge. State-of-the-art optimization techniques are for example 
sometimes benchmarked on their capability to find the global energy minimum of Lennard- Jones 
clusters.-"- Simulation studies have been devised, using enhanced sampling methods, for example 
aggregation- volume-bias Monte Carlo and umbrella sampling, to obtain the free energy and density 
of state of clusters or crystalline nuclei in liquids.-"— The recent boost of nanoscience increasingly 
demands a better and better ability to characterize, calculate, and understand theoretically the 
peculiar properties of matter at the nanoscale, a subject of broad technological, practical as well as 
conceptual implications. 

The traditional manner to approach cluster simulation is quite simple. As a first step one at- 
tempts to identify all the locally stable low energy structures, usually by some global optimization 
technique, such as simulated annealing, genetic algorithms, etc . The computational cost of these 
techniques is largely dominated by repeated potential energy evaluation, and one must carefully 
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choose a level of description of the mechanics of atomic coordinates (quantum or semiempirical or 
classical) and of the electronic degrees of freedom allowing at the same time an exhaustive search 
and a sufficiently accurate estimate of the energy. Sometimes a preliminary search with a less 
accurate empirical intermolecular potential can be followed by a subsequent refinement leading 
to a more accurate energy estimate and thus to the identification of the relevant structures with a 
higher level description. If one is interested in the properties of the system at a very low tempera- 
ture, what matters is the global energy minimum. When however temperature grows so that ksT is 
comparable to the energy gap Aqi above the global minimum, higher energy states become relevant 
and must be identified. It is traditional to consider a relatively small number of "relevant states", 
analogous to inherent states in glass physics, crudely consisting of locally stable atomic configu- 
rations or cluster geometries, that identify local energy minima in the space of atomic coordinates. 
A finite set of relevant states can be obtained for example by (real or conceptual) annealing of 
the system from a much larger number and variety of finite temperature instantaneous molecular 
configurations, cooled down to reach the closest local energy minimum at T = 0. At finite temper- 
ature, each relevant state s represents the center of a basin in configuration space. The logarithm 
of the volume of the relevant state basin volume in configuration space represents its entropy 5^. 
When the thermal occupancy of each relevant structure s is explicitly calculable, this is equivalent 
to specifying its free energy Fs = Es — TSs . A crucial question that needs to be addressed is, what 
contributes to the entropy of the relevant state, and how it can be properly calculated. The simplest 
and commonest approach is to consider free rotations and harmonic vibrations around the T = 
local energy minimum of the relevant state, and calculate the vibrational and rotational entropy 
involved, in the rigid-rotor-harmonic approximations.-^^ This is however not generally ade- 
quate, and in fact the rot- vibrational approximation may yield free energies of relevant structures 
in serious error. The full entropy of a state, measuring the size of its stability basin should include 
not just vibrations and rotations, but all possible coordinate transformations leading to essentially 
equivalent configurations. Proper determination of configurational entropy is a factor which crit- 
ically reflects the abundance of a relevant state in the cluster's thermal equilibrium. Although 
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these concepts are natural, and well rooted in classic cluster literature,-*^ they are otherwise not 
always properly implemented in everyday practice where configurational entropy, difficult to de- 
fine and evaluate, is often overlooked. A practical and accurate calculation approach including 
automatically all sources of entropy, especially for the more complex molecular clusters, is highly 
desirable. 

In this paper we introduce and demonstrate a procedure leading to a computationally conve- 
nient and conceptually correct extraction of the free energy of the relevant states of a cluster at finite 
temperature, including all non-electronic entropy sources contributions. In the simple case we use 
for our demonstration, a small water cluster, the entropy is found to differ drastically between one 
relevant state and another. This working case also shows explicitly that the neglect of configura- 
tional contributions leads to large errors, totally unacceptable practically no less than conceptually. 
Our method simultaneously addresses both steps normally followed in traditional cluster study 
techniques: a) the search for the relevant structures, and b) the study of their finite- temperature 
properties. The approach, based on a recently developed enhanced sampling technique, replica 
exchange metadynamics,— allows the computation of the true free energy of each relevant struc- 
ture, directly sampling its basin and measuring its occupancy (overall abundance) in full thermal 
equilibrium, thus including automatically all entropy sources. 

For our demonstration purposes we chose the water nonamer (1120)9, a very simple dielec- 
tric cluster which nonetheless displays a sufficient complexity and interesting physics in its rele- 
vant structure spectrum. Although not directly addressed in this study, well characterized spec- 
troscopic properties of this particular cluster are available.— We studied the cluster properties 
covering a temperature range of potential interest for atmospheric science (100-250 K), and de- 
scribed the forces between water molecules with a simple but realistic and widely tested empirical 
force scheme, the so-called TIP4P potential.— Since the electronic shells of this system are closed 
and widely gapped, the use of such effective interatomic forces is not unreasonable, although 
clearly not definitive. The cluster thermal evolution is simulated by replica exchange metady- 
namics (RE-META),^"^ an approach which requires the simultaneous running of several molecular 
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dynamics simulations at different temperatures. Each simulation ("replica") is biased by an ad- 
ditional history-dependent "metadynamics potential".— This artificial bias energy term acts on an 
appropriately chosen small set of collective variables (CV), global or local order parameter-like 
quantities representing some physically motivated function of molecular coordinates. At fixed 
time intervals, the replicas are forced to attempt an exchange of their coordinate configurations. 
The exchange move is accepted or rejected according to the usual Metropolis criterion,—^ in- 
cluding the change of the history dependent CV-based bias potential between the two replicas. 
The configurations visited in the RE-META simulation correspond to visiting the free energy land- 
scape of the cluster - as a function of the CV variables - in parallel at all temperatures. Combining 
replica exchange and metadynamics leads to an overall kinetics which, even at low temperatures, 
rapidly becomes diffusive in collective variable space, resulting in an extremely efficient configura- 
tion sampling.— In RE-META, replica exchange is improved because metadynamics explores high 
free energy configurations within each replica, ultimately leading to faster decorrelation. Metady- 
namics is in turn strengthened, since replica exchange improves the sampling of the degrees of 
freedom that are not explicitly biased - thus mitigating one major problem of metadynamics.— 
The water nonamer example confirms that the free energies of the lowest relevant states are indeed 
heavily affected by configurational entropies. In particular, the entropy of the first excited state is 
higher than that of the ground state, which causes for increasing temperature an outright crossing 
of the two free energies and populations, a crossing which does not occur in the vibrational and 
rotational approximation. 



The (H20)9 cluster: RE-META calculations 

For the water nonamer (H20)9, we carried out RE-META simulations with eleven replicas at the 
following temperatures: 25.8, 35.8, 47.8, 62.9, 82.4, 108.3, 136.4, 170.0, 205.5, 248.7 and 298.4 
K. Since TIP4P bulk ice has a melting temperature of 232 K at 1 bar,— the higher temperatures 
correspond to a fast diffusing liquid-like regime, the lower ones to a slow diffusing solid-like one. 
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Our level of description is entirely classical, and no quantum effects are included, either nuclear 
or electonic. Because therefore our treatment ignores all quantum freezing and zero point motion 
effects, the simulation results lose their strict validity at very low and zero temperatures, a limit 
where they should not be extrapolated. The molecular dynamics simulation was performed in 

o 3 

the NVT ensemble and a cell volume of 18850 A , a volume where the probability of observing 
an isolated water vapour molecule is extremely small. A Nose-Hoover chain thermostat was used 
to control the temperature, the integration time step was set to 0.2 fs, and configurations were 
stored for analysis every ps. Each replica was evolved for a total simulation time of 50 ns and 
independently biased by a metadynamics potential acting on a single collective variable S counting 
approximately the number of hydrogen bonds inside the cluster: 



1 - inj /roY 
ij 1 - injf 

where the sum runs over proton-oxygen pair bonds belonging to different molecules and ro = 2 A. 
The metadynamics potential has the form— 



< t 



where s{t) — S{r{t)), 5s — 0.15 and w — 0.004 meV. As customary in metadynamics, the 
collective variable defined in ?? is chosen on purely physical grounds, and is neither unique 
nor systematically definable in a general case. However, while this arbitrariness could be seen 
as endangering the reliability of metadynamics, the replica exchange strategy does, as shown in 
Ref., effectively cure the problem. In our calculation the exchange moves between configurations 
taken at different temperatures were attempted every 2 ps, accepting the moves with a probability 
min(l,exp(— A)) with 
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where ij is the index of the rephca, V is the sum of the potential energy and the bias, r are the 
coordinates and jS is the inverse temperature, which is different in each rephca.— We found that 
after a typical equilibration time Tf ^ 4 ns the history-dependent potential started growing evenly 
in CV space - meaning independent of 5 - in all replicas, indicating that -VciS^t) will eventually, 
for large t, provide a good estimate of the free energy as a function of S.— After the equilibration 
time Tf we could therefore stop changing the history-dependent potential Vg and start collecting 
the statistics of occurrence of properly defined relevant states, with the scope of computing their 
respective probability, and from that their free energies and entropies in thermal equilibrium. 

In order to analyze the results of the procedure it is necessary to assign during the MD simula- 
tion each molecular coordinate configuration, to the closest reference structure - its relevant state. 
If all the possible local minima of the potential energy surface were known, one could classify the 
configurations according to their distance from the reference structures, measured, for instance, 
by the appropriate root mean square deviation of coordinates. This however would require a pre- 
liminary search of all the minima. More importantly, in order to assign a configuration to the 
correct relevant state, one should consider all permutations between identical atoms, a procedure 
hugely cumbersome and computationally expensive. We apply here a more conservative procedure 
that uses ideas from the basin hopping technique.- Starting from all the coordinate configurations 
generated by MD, we carried out a static potential energy minimization by conjugate gradients 
evolving atomic coordinates until the force on each atom was less than 0.0001 eV/A. If a reference 
structure with the same potential energy is already available, the instantaneous MD configuration is 
assigned to that structure. Otherwise, a new reference structure is introduced in the pool. Here two 
quenched structures are assumed to be the same if their energy difference is smaller than 0.0001 
eV. This procedure is still computationally intensive, but allows an unambiguous assignment of all 
MD frames and is insensitive to permutation of the atoms. 



[figure] [!][]! shows the number of independent structures that are explored as a function of 



simulation time at different temperatures for a normal replica exchange run, which we carried out 
in parallel for comparison, and for RE-META. The number of structures visited is much larger 
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in the latter, a beneficial result of the metadynamic bias technique, especially important at low 
temperature. This demostrates how RE-META can be successfully used to carry out the first task 
normally undertaken in cluster studies namely the search for all the low energy relevant structures. 
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Figure 1 : The number of relevant states found as a function of simulation time at different tem- 
peratures with simple replica exchange (RE) and with RE-META. Note the great improvement 
brought by addition of the history dependent metadynamics potential ( ??), especially important at 
low temperature. 



[figure] [2] []2| -a shows the molecular structure of the water nonamer lowest energy relevant 
structures. The five lowest energy structures share the same oxygen arrangement, and differ only 
by the position of the protons. Eight oxygens form a regular hexahedron and an additional oxygen 
is attached to one of the edges. The global minimum is identical to that previously identified in 
ref In all the nine relevant structures 13 H-bonds are formed. 



As seen in [figure] [2] []2 the lowest energy structure has a mirror symmetry about the plane 



shown as a black line in [[figure] [2] []2| -b, while the same mirror plane gives rise to two distinguish- 
able degenerate structures for the second lowest structure. This, as we will see in the following, 
significantly influences their respective configurational entropies. 

We now come to our main point. Besides providing a tool for the quick exploration of cluster 
structures, RE-META crucially yields an accurate estimate of the probability to observe each rele- 
vant structure 5 as a function of temperature. Following Ref, we first compute the average bias 
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Figure 2: (a)Lowest energy relevant states of the water nonamer. (b) Front and back viewed 
configuration for state (lowest energy state). Black line: symmetry plane, (c) Front and back 
viewed configurations for state 1 (second lowest state). For a fixed oxygen skeleton state 1 has 
twice larger multiplicity than state 0, which has a mirror symmetry. Thus the configurational 
entropy of state 1 may be expected to exceed that of state by approximately In 2. 



potential 



VG{S) = l/{Tsim-TF) dtVoiS^t), (4) 



where Tf is the equilibration time of Vg, 4 ns in this work and T^/^ is the total simulation time. 
Denoting by r/ all the configurations generated by RE-META at temperature T, we estimate the 
equilibrium probability to observe relevant structure a as 

Y,oxp{VG{Siri))/kBT) 
p"(T) = — — (5) 

where the notation / G a means that the sum runs over the configurations that are assigned to struc- 
ture a as explained above. In order to assess the reliability of this approach, we also computed the 
probability pa in normal replica exchange where the probability is simply proportional to the num- 
ber of times the structure a is observed. As shown in |[figure] [3] []3f a,b, the probabilities estimated 



in the two manners are consistent. The small remaining discrepancy between RE and RE-META 
is due to statistical uncertainties. We remark that the errors of methods based on statistical sam- 
pling are inevitably large at low T, as the relevant transitions are observed relatively rarely even if 
a powerful enhanced sampling approach is used. 
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Figure 3: (a) occupancy probability for the lowest energy relevant state ("0") versus tempera- 
ture, with RE-META, normal RE, rigid-rotor-harmonic approximation ("Rot-Vib") and rigid-rotor- 
harmonic approximation with an additional the In 2 factor in the entropy of state 1 - state 8 ("Stan- 
dard"), (b) occupancy probability for the second lowest energy relevant state ("1"). The statistical 
error bars in RE and RE-META represent the standard deviations estimated by separately com- 
paring averages in three successive portions of the simulations. Note the substantial agreement 
between RE and RE-META and the large deviation of the rot-vib approximation, showing its clear 
inadequacy. Note also that the occupancies of states computed with RE and RE-META are fully 
consistent with those computed using the standard approach at low temperature, and how the prob- 
ability of state 1 overtakes that of state around 150 K, as a result of its larger configuration 
entropy. 
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Comparison between RE-META and rigid-rotor-harmonic vi- 
bration approximations 



We calculate and compare in [[figure] [3] []3f a,b the exact occupancy probabilities of the two low- 



est relevant structures with the approximate probabilities, now estimated using the rigid-rotor- 
harmonic vibration approximation to the cluster vibrational spectrum of each state: 

exp(-F«A^r) 
^ ^ ' Lpcxpi-FP/kBT) 

where the rigid-rotor-harmonic free energies are calculated as 

F«=£«-^rin^«,_,,, (7) 

where is the energy and 

^rot-vib ^ ^rot + ^vib- 

^- = n^- (10) 

where /-^ are the principal moments of inertia, is the rotational symmetry number and 0)^ is 
the vibrational frequency of mode m obtained by diagonalizing the Hessian matrix of the potential 
in structure a at r=0. To keep the calculation simple, the summation in the denominator of ?? 
runs over the nine lowest energy states. We have checked that if 30 structures are included in the 
summation, the probability to observe all states between the 9th and 30th is only 0.2% at 108K and 
13% at the highest temperature 250K. It should be noticed that the rotational symmetry number 
of the lowest energy relevant states of the water nonamer are all equal to 1 and the moments of 
inertia are similar. So, for this system, the rotational entropy varies much less than the vibrational 
entropy. For example the rotational entropy difference between state and state 1 is only 7 % of 
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the vibrational entropy difference. As is seen in |[figure][3][]3[ the true occupancy probabihties 
from RE-META and those predicted by the rigid-rotor-harmonic vibration approximation ?? differ 
very significantly. For instance, at T = 200 K, ?? still predicts a significant population for state 
and state 1, while RE-META shows that the population of state has become negligible, while 
the population of state 1 is only about 10 %. Moreover, the rot-vib approximation ?? implies that 
state is always more populated than state 1, whereas according to RE-META the most populated 
state at 150 K is state 1 and not state 0. 

We finally calculated the entropy difference between states and 1 

^E-^F 

A5= (11) 



where AF and lS.E are the free energy and the internal energy differences, [[figure] [4] []4| shows 
the entropy difference estimated with RE-META, RE, and rot-vib approximation. This entropy 
difference is finite and of order 1 even at the lower temperatures signalling an initial configurational 
degeneracy of state 1, and grows further upon heating. The rot-vib approximation on the other 
hand misses both effects. Visual inspection of the structures assigned to relevant states and 1 
(data not shown) reveals that up to T=200 K the oxygen skeleton remains largely unchanged. Thus 
the main source of configurational entropy is proton disorder. 

The state has a mirror symmetry, while state 1 does not. The broken symmetry implies that 
state 1 is two-fold degenerate. Thus, the relative configurational entropy of state 1 is roughly In 2, 
which, at low temperature, approximately agrees with the relatively large occupancy of state 1. 
Adding In 2 to the rot-vib entropy for states (but not for state 0) leads to a positive "Standard" 



entropy difference (black curve in [[figure] [4] []4) ), closer to the simulation results at small T . At 



high temperature, anharmonic contributions further favor state 1, making the entropy difference 



estimated even with the standard approach unreliable. Consistently, as shown in [[figure] [3] [] 3 [ the 
occupancies estimated with the standard approach are compatible with those estimated with RE- 
META at low r, but important deviations appear at T > 100^. The comparison of the occupancies 
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presented in ?? and, even more, of the entropy differences of [[figure] [4] []4[ are the central result 



of this work. 
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Figure 4: Entropy difference between state 1 and state as a function of T. Note the initially larger 
configurational entropy of state 1, further growing with temperature, due to anharmonicity. Both 
effects are missed by the rot-vib approximation, whereas the addition of ln(2) (the "Standard" 
approximation) provides a reliable result only well below 100 K. 



Discussion and conclusions 

The cluster description obtained with RE-META can naturally be used to predict measurable quan- 
tities. As an example, we computed the average dipole of the water nonamer as a function of T. In 
[figure] [5][]5[ the accurate RE-META computed dipole is compared with the rot-vib approxima- 



tion: D — Y.a^^P^ (^) where is the dipole in the reference structure a and the probabilities 
Pa are estimated by ??. The exact and the approximate results differ very significantly, underlining 
the devastating effects of ignoring the configurational and anharmonic entropies, especially at high 
T. 

Beyond the deliberately simple case study (H20)9 presented here, RE-META method is appli- 
cable to larger water clusters, and of course to different systems, once the collective variables are 
properly selected. The workload and complexity grows considerably with size - not surprisingly 
in view of the exponential growth of phase space - but the effectiveness is not impaired. As an 
example we considered the larger water cluster (H20)^4, hundreds of times more complex than 
the nonamer with a 10 ns simulation, [[figure] [6] []6| shows the number of independent structures 



explored as a function of simulation temperatures for RE and RE-META, underlining the strong 
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Figure 5: Thermal averaged dipole moment of a water nonamer, computed directly from the RE- 
META and RE trajectories and using the probabilities ( ??) estimated within the "Rot-Vib" and 
"Standard" approximation. 

superiority of the latter, in particular its ability to unearth relevant structures that would otherwise 
be simply ignored. 
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Figure 6: The number of relevant states of (H20)^4 found as a function of temperature with simple 
replica exchange and with RE-META. Improvement brought by RE-META method is also notable. 
Insets: (1) Lowest relevant state. (2) Second lowest relevant state. 

In summary, we have shown how the full configurational complexity of a small cluster can 
be directly, realistically and economically accessed by a suitable sampling, combining replica ex- 
change with metadynamics techniques (RE-META). This proves extremely effective in searching 
the cluster's relevant states, and in measuring their thermal population as a function of T. In the 
water nonamer, these populations are shown to differ largely from those routinely predicted by 
rot-vib approximation. The crucial factor which influences the populations is the configurational 
entropy of the relevant states, an entropy which is generally inaccessible but which in RE-META 
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is naturally and directly measured and understood. The method holds a good promise of applica- 
bility in a great variety of relevant systems, provided a set of good, physically motivated collective 
variables can be identified. 
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